#include "FEMSpace.h"
#include "BoundaryFunction.h"
#include <math.h>
#include <ctime>
#include "Error.h"
#include <fstream>
#define pi 4.0 * atan(1.0)

double u(double *p)
{
    return sin(pi*p[0])*sin(pi*p[1]);
}
double u_x(double *p)
{
    return pi*cos(pi*p[0])*sin(pi*p[1]);
}
double u_y(double *p)
{
    return pi*sin(pi*p[0])*cos(pi*p[1]);
}
double f(double *p)
{
    return 2*pi*pi*sin(p[0])*sin(p[1]);
    //return 2 * pi*  pi* cos(pi * p[0]) * cos(pi * p[1]);
}

double bc(double *p)
{
    return sin(pi*p[0])*sin(pi*p[1]);
    //return cos(pi * p[0]) * cos(pi * p[1]);
}
int main()
{
    
    //g++ -o main test_possion.cpp -std=c++11 -I /usr/include/eigen3/
    RectangleDomain* r = new RectangleDomain({{0,0},{1,0},{1,1},{0,1}});
    std::vector<Boundary<2> > B = r->boundary();
    int n,order;
    std::string elementtype;
    std::cin >> n >> elementtype >> order;
    Mesh<2> *m;
    Element<2> *e;
    if(elementtype == "P"&& order ==1)
    {
	m = new P1Mesh(r, {n, n});
	e = new Triangular_1_Element();
    }
    else if(elementtype == "P"&& order ==2)
    {
	m = new P2Mesh(r, {n, n});
	e = new Triangular_2_Element();
    }
    else if(elementtype == "P"&& order ==3)
    {
	m = new P3Mesh(r, {n, n});
	e = new Triangular_3_Element();
    }
    else if(elementtype == "P"&& order ==4)
    {
	m = new P4Mesh(r, {n, n});
	e = new Triangular_4_Element();
    }
    else if(elementtype == "Q"&& order ==1)
    {
	m = new Q1Mesh(r, {n, n});
	e = new Quadrilateral_1_Element();
    }
     else if(elementtype == "Q"&& order ==2)
    {
	m = new Q2Mesh(r, {n, n});
	e = new Quadrilateral_2_Element();
    }
     else if(elementtype == "Q"&& order ==3)
    {
	m = new Q3Mesh(r, {n, n});
	e = new Quadrilateral_3_Element();
    }
     else if(elementtype == "Q"&& order ==4)
    {
	m = new Q4Mesh(r, {n, n});
	e = new Quadrilateral_4_Element();
    }
    PoissonEquation* equ = new PoissonEquation();
    equ->SetRightHandsTermFunction(f);
    BoundaryFunction<2>* dir = new Dirichlet<2>(bc,B);
    BoundaryCondition<2> bfc(dir);
    Poisson_2D *poissonproblem = new Poisson_2D(m, e, equ, bfc);
    poissonproblem->AssembleStiffMatrix();
    poissonproblem->AssembleRightHandsTerm();
    poissonproblem->DealWithBoubdaryCondition();
    poissonproblem->Solve();
    
    Eigen::VectorXd solution = poissonproblem->solution();
    //std::cout << solution<<std::endl;
    std::cout << n << "L2_error: " <<  poissonproblem->Error_L2Norm(u) << std::endl;
    std::cout << n << "H1_error: " << poissonproblem->Error_H1Norm(u,u_x,u_y)<<std::endl;
     std::ofstream output("./output/result2.txt");
     for(int j=0; j <solution.size();j++)
     {
    	 Dofs<2> dof = m->DofsofIndex(j);
    	 output << dof[0] << " "<<dof[1] << " " << solution[j] << std::endl;
    }
    output.close();
    // std::vector<int> a = m->NodeofEle(2);
    // for(int i= 0; i < a.size();i++)
    // {
    // 	std::cout << a[i] <<std::endl;
    // }
    return 0;
}
